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Abstract. We study structural transitions in a system of interacting particles 
arranged as a crystalline bilayer, as a function of the density p and the distance 
d between the layers. As d is decreased a sequence of transitions involving 
triangular, rhombic, square and centered rectangular lattices is observed. The 
sequence of phases and the order of transitions depends on the nature of 
interactions. 



1. Introduction 

Structural transitions in solids may be caused by various external parameters such as 
temperature, pressure, stress, electrical and magnetic fields [1, 2jetc. Confinement and 
dimensional reduction can also lead to structural transitions especially in soft solids 
like colloidal suspensions. Colloidal solids are especially suited to modification and 
manipulation using a variety of means such as structural confinement pj , laser-induced 
phase transitions 0, shear 0, static [6] and dynamic [7] external fields. 

Confined colloids kept in a thin wedge geometry of two optically-flat quartz glass 
plates, exhibit a sequence of structural transitions : nA (n + l)n (n + l)A with 
increasing wedge height, where n is the number of layers and A and □ corresponds 
to layers of triangular (p6) and square(p4m) symmetry respectively [8^, '9 . The full 
equilibrium phase diagram of such a system has been studied analytically as well 
as using extensive computer simulations [101 HI]- The transitions are usually first 
order, though continuous transitions via a layer buckling mechanism has also 
been predicted and observed. 

In this paper, we explore another way in which structural transitions may be 
induced in a colloidal solid. Consider a crystalline bilayer separated by a distance 
d between the layers. Each of these layers is held in place by individual trapping 
potentials, set up, for example using laser tweezers The strong trapping potential 
ensures that out of layer fluctuations are typically unimportant. We investigate the 
stability of the bilayer as the distance d is decreased. We show that d behaves as a 
controlling parameter and induces a rich sequence of transitions involving a variety of 
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two-dimensional Bravais lattices. The exact sequence of transitions crucially depends 
on the nature of interactions. In this paper we study two kinds of model solids (a) the 
generic Lennard Jones (LJ) [T4] solid and (b) the soft Gaussian core model (GCM) 
appropriate for suspensions of globular polymers. Our main results are as follows. For 
the LJ system, we obtain at temperature T = two independent triangular (TRN) 
crystalline layers for large d. As d is reduced, the system undergoes a first order 
transition to a staggered square (SQR) solid. As d is further reduced, this square 
solid becomes, first, a centered rectangular (CR) and finally again a triangular solid 
as d is decreases to zero and the layers merge. The final two transitions are continuous. 

In contrast for the GCM, all the transitions are continuous. As d reduces, the 
TRN solid transforms to a SQR solid continuously through a sequences of rhombic 
(RMB) lattices. The SQR solid subsequently transforms back to the TRN solid for 
small values of d again continuously but this time it uses a sequence of CR lattices with 
intermediate aspect ratios. The progression of phases seen correspond roughly with 
those seen in the extensive literature on classical interacting bilayer Wigner crystals 
[16], though the sequence of phases and the nature of transitions are different. 

The rest of the paper is organized as below. In the next section, we describe 
the bilayer system in detail, introducing the order parameters for the transition and 
state the interatomic potentials used. In Section 3,we give the results for the zero 
temperature energy minimization. This is followed in Section 4 by a full normal mode 
analysis investigating the stability of the ground states obtained in Section 3 and 
the nature of the transition. In Section 5, we present results of finite temperature 
Monte Carlo simulations. We discuss our results and their implications and conclude 
in Section 6. 

2. Model System 



Figure 1. A schematic diagram of the model bilayer solid explaining the 
structural parameters. The crystal structure on each of the layers, denoted by 
filled and open circles, are identical but staggered by the amount A along the 
diagonal. The lattice parameetrs are a and b and the apex angle is cf). The half 
filled circle is the projection of a lattice point (filled circle) on the top layer to the 
bottom one. 



d 
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Consider a system of 2N particles arranged as two parallel two-dimesional 
crystalline layers of N particles each (see Fig.([T])). The crystal structure of each 
of these layers may be assumed to be a general two-dimensional oblique {p2) lattice. 
Each particle interacts with all other particles via isotropic and pair-wise interacting 
potentials. The position vector for the i-th particle of this lattice can thus be expressed 
as 



where m, n = 0, 1, 2, • • •, a and b are two in-plane basis vectors, (j) is the angle between 
these two basis vectors, d is the interlayer separation and A is a shift between the 
center of masses of two layers (Fig.([T])). Therefore to specify our model completely, 
knowledge of these five variables, a, 6, (/), A is sufficient. The last variable A ensures 
that for small values of the particles from different layers do not overlap. The 
particles are not allowed to fluctuate out of the layers and particle exchange between 
the layers is prohibited. Each layer is therefore considered to be strongly confined in 
the Z-direction while they are allowed complete freedom in the X,Y- plane. 

Let us first imagine the possible physical scenario as the layers are brought close 
to each other starting from a large inter-layer separation. When the layers are well 
apart they exist as two independent mono-layers and they show TRN symmetry, 
the minimum energy configuration for the two-dimensional crystalline system for the 
interaction potentials considered by us. As the inter-layer separation, between 
these two layers starts decreasing, the system passes through a series of structural 
transitions which may involve RMB and SQR phases. For even smaller values of d, 
two layers start merging into one and the TRN symmetry is regained. Transformation 
between a SRQ and a TRN phase may be accomplished, in general, by either (i) shear 
i.e. change in the angle between two in-plane basis vectors producing an intermediate 
RMB structure or (ii) change in the aspect ratio (b/a) which produces a CR lattice. 
Both the RMB and the CR lattices being less symmetric have TRN and SQR phases 
as limiting cases. 

It is therefore clear that we need to introduce two order parameters [17 in order to 
describe completely the phase transitions of our model system. First, the bond angle 
order parameter ip = coscj) which is when system takes SQR symmetry {(j) = 90°) 
and non-zero otherwise. The second order parameter ^ is related to the aspect ratio 
b/a as C = {b/a - 1)/(V3 - 1) which varies from in the SQR to a non-zero value 
in the CR phase. Note that the highly symmetric TRN phase is described both by 
(V^ = 0.5, ^ = 0) and {ip = 0.0, ^ = 1). Finally, if Cij^i^j = x^y is the two dimensional 
strain tensor, the shear, 63 = Cxy^ and deviatoric, 62 = Cxx — ^yy^ strains are related 
to and ^ as. 



We have studied phase transitions in the bilayer system for two different model 
potentials. The Lennard Jones potential : 



fi = (m + A)a + (n + X)b + dz 

= ((m + X)a + (n + X)b cos (j))x + (n + X)b sin (j)y + dz 



(1) 



tan0 = 63/(1 — ei) 

2(l + e2)+ei 
^ 2(1 - £2) + ei 



(2) 




(3) 
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has been used extensively in the past as a generic model which includes both long 
range attractive and short range repulsive interactions. In Eq. ([3]), r^j = lu — 
the distance between the i and j-th particle. An intrinsic length scale rmin = 2^/^<j 
corresponding to the minimum of Vlj niay be associated with this potential. The 
nearest neighbour distance between particles is close to this value throughout. We use 
reduced units for LJ potential throughout the paper defining lengths in units of a and 
energy in units of e. It follows that the densities are in units of and temperatures 
in units of e/Zc^, where /cb is the Boltzmann factor. 

On the other hand, Gaussian core potential has been used specifically to model 
soft solids. 

U= Y^VocMinj) 

= J2Voexp[-{rij/l)% Vo,l>0 (4) 

Vq and / acts as the energy and length-scales for this potential and r^j is same as it was 
for the previous equation. The GCM is interesting because, firstly, the potential is soft 
and is known to display behaviour similar to that of real polymeric solids. Also many 
of its properties are well known and tested especially because at very low densities it 
reduces to the hard disk model which is widely used for modelling colloidal solids. Due 
to the purely repulsive nature of this potential, it does not have any preferred nearest 
neighbour distance which is determined in this case by the density. This potential 
also possesses an interesting duality property [18j such that high (p>) and low (p<) 
density properties are related to each other by 

P>P< = TT-^ (5) 

where D is the dimensionality. It is therefore sufficient to confine our studies in the 
range < p < p/.p. where the fixed point density pj.p. = tt"^ 0.32 in two-dimensions. 

3. Zero Temperature Calculation 

In this section, we determine the ground states of our system in the space of the order 
parameters, (j) and ^, keeping the layer separation d and the density p as external 
parameters. Below we discuss our results for the LJ solid and the GCM one after the 
other. 

3.1. Lennard- Jones potential 

We begin with 2N particles divided into two layers each arranged in a TRN lattice 
of N particles. For fixed layer separation and the lattice parameter a set by the 
density p, we minimize the total energy U with respect to A and b. The resulting 
order parameter ^ remains constant at 1 - the value appropriate for the TRN lattice. 
For the lattice sums needed to calculate the minimized energy we have set a cut-off 
radius, Tc = 2.5cr. 

In Fig. ([2^) we have plotted minimized energy Umin as a function of the order 
parameter ip for various values of d. For large inter-layer separation, bond angle order 
parameter at the minimized energy (Umin) shows a non-zero value, i/j = coscj) = 0.5, 
which indicates TRN symmetry {(f) = 60°). The minimum in A is very shallow 
indicating two triangular lattices prefer to remain independent. As d is decreased 
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Figure 2. (a) Plots of the total energy as a function of -0 at p = 1.2 and (top 
to bottom) d = 0.66, 0.68, O.TOand 0.72 showing the SQR to TRN transition for 
the LJ system, (inset) The energies of the SQR (filled square) and TRN (filled 
triangles) plotted as a function of d are seen to cross with a change in slope as 
expected of a first order transition. The fact that this transition is first order has 
also been verified by a normal mode analysis, (b) Plot of the order parameter 
^ as a function of d at the same density showing the two continuous transitions 
SQR ^ CR ^ TRN as d ^ 0. 



two layers fall into registry and for further decrease of d a second minimum at = 
develops corresponding to the SQR phase {<p = 90°). A plot of the energy of the SQR 
and TRN phase (Fig. ([2^) inset) shows a first order transition at d = 0.684. This fact 
is confirmed by a normal-mode analysis presented in the next section. 

At even smaller values of the square solid begins to deform continuously by 
changing the aspect ratio (or order parameter ^) at fixed ^. In Fig. ([2)3) we have 
plotted the minimized value of ^ vs. d. From the plot it is obvious that there are 
two distinct continuous transitions taking the solid from SQR at ^ = to an eventual 
TRN phase at ^ = 1 via an intermediate CR lattice of ^ = 0.5. 

3.2. Gaussian Core Model 

Due to the purely repulsive nature of the Gaussian core potential, it does not have any 
preferred nearest neighbour distance and at zero stress the solid would disintegrate. 
We have carried out all minimizations for a system with fixed reduced density, p = 0.2. 
We have plotted the minimized energy Umin as a function of the bond angle order 
parameter cj) for various values of d. In Fig.([3ti) in contrast to the Lennard- Jones solid 
we now obtain continuous transitions from TRN to SQR through a set of RMB phases 
and back to TRN at d = through a set of CR phases. 

3.3. T=0 Phase Diagram 

Our results concerning the various ground states and structural transitions in both 
the LJ and GCM systems in the p - d plane has been shown in the zero temperature 
phase diagrams, Fig.(|4^) and (b) respectively. Note that for both the systems the 




Figure 3. (a) Plot of the total energy in the GCM for p = 0.2 and (top to 
bottom) d = 1.4, 1.5, 1.6, 1.7 and 1.9 as a function of ip. Note the appearance of 
the minimum at -0 = 0.5 corresponding to the TRN lattice as d increases, (b) 
The continuous transition from the SQR to the TRN lattice as ci ^ 0. 




Figure 4. Zero temperature phase diagrams in the p — d plane for the LJ (a) 
and GCM (b) systems. The various phases are marked. First order transitions 
are shown by solid and continuous transitions by dashed lines. 



triangular phase is stable at all p for both very large d and (i = where the system 
becomes effectively two dimensional. 

4. Normal Mode Analysis 

To further elucidate the nature of the structural transitions in the two systems, we 
have undertaken a normal mode analysis [19j of the TRN and SQR solids obtained 
for each of the two interactions. 
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Let u{fi) be the displacement of the i-th particle from its equilibrium position r^. 
Within harmmonic approximation, now the potential can be written as 

Uharm ^^u{fi)V{fi - fj)u{fj) (6) 

V(ri-fj) = D^^{ri-rj) 

= ^r^rJ ^ y^ivifi - fk) - V^yifi - fj) (7) 
k 

where V^zy(r) = d^V j dr ^dry. We have 2N equations of motion, one for each of the 
three components of the N particles, since we have already restricted fluctuation in 
the z-direction. 

^(^0 = - - r,)u{f,) (8) 

3 

We seek solutions to the equations of motion in the form of simple plane waves : 

u{f^t)=e exp[i{k.f — ujt)] (9) 

Here e is the polarization vector of the normal mode. The Born-vonKarman periodic 
boundary conditions restricts the wave vector ^ to a single primitive cell of the 
reciprocal lattice vector, which is normally identified with the first Brillouin zone. 

Substituting Eq.([8]) into Eq.Q we find a solution of the three-dimensional eigen 
value problem : 

u;^e = V{k)e (10) 
Here V{k)^ the dynamical matrix, is given by 

V{k) = ^V{n)e-'^^-^'' (11) 

i 

Two solutions to Eq.(|9]) for each of the N allowed values of k give us 2N normal 
modes. The reciprocal lattices for both SQR and TRN are known to have the same 
symmetry of the real space lattice. Exploiting this property, one finds the values of 
LO only for those /c- values along the lines connecting the high symmetry points of the 
first Brillouin zone (Fig. ([5^)) thereby obtaining the dispersion curve a; vs \k\ and mode 
structure of a given lattice. For any stable equilibrium structure co should always be 
non- negative definite. 



4.1. The Lennard Jones potential 

To show that the SQR ^ TRN transition at large d is indeed first order, we have 
obtained the dispersion curves for the metastable SQR phase for a value of d slightly 
larger than the critical dc = 0.684 for the chosen p = 1.2. This is shown in Fig. ([5)3). 
This indicates that the transition is first order with the possibility of co-existence. 
When d is decreased further, at the minimized energy Umin shows minima at = 
but the normal mode analysis shows us that the SQR structure can not be stable at 
this value of d {Fig.^)). Actually, the two layers start merging into one, by changing 
the aspect ratio away from b/a = 1. 
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Figure 5. Results for phonon dispersion curves : (a) diagram showing the high 
symmetry points in the SQR and TRN (reciprocal) lattices, (b) Phonon dispersion 
for the stable SQR lattice at p = 1.2 and d = dc = 0.684 i.e. the value of d at SQR 
to TRN transition. Note that the dispersion curve shows that the SQR phase is 
locally stable pointing to a first order transition, (c) Phonon dispersion of the 
SQR solid at p = 1.2 and d = OA showing an instability in the transverse acoustic 
branch. At this value of d the CR solid is stable, (d) Dispersion curve for the 
TRN solid in the GCM for p = 0.2 and d = 1.6 showing the appearance of a non 
zero k vector soft mode. 



4.2. The Gaussian Core Model 

For the Gaussian core model the scenario is quite different. For intermediate values 
of (i, the SQR structure is seen to be unstable and the mode structure for the 
TRN structure exhibits mode softening, Fig.([5]i)), such that u for /c 7^ 0. 
Examination of the deformation corresponding to this k shows that the SQR lattice 
becomes unstable to shear deformation at the zone boundary. This mode softening 
therefore establishes the transition from TRN to SQR transition to be continuous for 
the Gaussian core model. 

The transition from SQR back to TRN at low d is always continuous both for LJ 
and GCM as verified by our normal mode analysis. 

5. Finite Temperature results: 

We close our discussion on structural transitions in a bilayer crystal by briefiy 
mentioning some of our results for the LJ case using Monte Carlo simulations [20j. A 
detailed calculation of the phase diagram of the bilayer GCM in the temperature, p 
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and d space using both Monte Carlo and classical mean field density functional theory, 
which is known to yeild particularly good results for this system, will be published 
elsewhere. 
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Figure 6. Configurations from our Monte Carlo simulation in the LJ system for 
T = 1.0. (a) and (b) show the result of increasing d from (i = 0.6to(i = 0.8at 
fixed p = 1.20. Note the first order SQR to TRN transition, (c) and (d) show the 
corresponding result for pressure induced SQR to TRN transition as p is increased 
from 1.3 to 1.4 at fixed d = 0.6. Our results are consistent with those seen in 
Fig-dSi). 

The simulation is done using the usual Metropolis algorithm keeping total 
number(A^) of particles, volume(V) and temperature (T) fixed. Periodic boundary 
conditions have been used for all directions except in the direction of stacking of the 
layers. For our purpose, we use a system of 512 number of particles, temperature is 
fixed at 1.0 and the volume is determined by the density, p = 1.2, used for simulation. 
Starting from a large layer seperation {d 0), we have observed the system change 
its symmetry from TRN to SQR with decreasing d. For large values of d {d > 0.8), 
when the two layers are well seperated, each layer behaves independently of each 
other exhibiting the expected TRN structure. As d is decreased, the two layers start 
interacting with each other. As a result, part of the system starts to transform into a 
SQR. Through this phase co-existence, the whole system transforms to a state where 
each layer shows SQR structure at d = 0.6. 

We have also explored the configurations of our model system in the temperature- 
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density plane keeping the inter-layer separation fixed. For d = 0.6, where we have 
already seen the SQR structure to be the ground state, we increase p to 1.3. We again 
encounter a first order boundary and the system equilibrates to the SQR structure. 

The values of the critical d and p compare favourably with the T = phase 
diagram shown in Fig. (|4^). Similar scans at other temperature values over a large 
range of p and d have confirmed that it is not possible to induce a structural transition 
in this model by changing T. The phase boundaries shown in Fig.(|4j) therefore extend 
vertically upto a melting temperature Tm{p^d). This behavior is identical to that 
seen in the classical Wigner crystal and is probably an universal feature of any such 
crystalline bilayer. 

6. Summary and Conclusion 

In this paper we have studied structural transitions in a bilayer crystal. We have shown 
that the system has a rich phase diagram and shows a number of phase transitons. 
The identity of the phases and the nature of the transitions depend on the interaction 
potential. On the other hand, some features of these transitions e.g. the overall 
topology of the phase diagram seems to be similar and independent of the details of 
the interaction. 

We believe that it may be easy to realize this soft matter system experimentally 
and study many of its interesting equilibrium and dynamic characteristics. For 
example, critical properties of the continuous structural transitions, especially, near 
the melting line and a detailed study of finite size scaling and crossover in these systems 
may be illuminating [21j. We are also particularly interested in the dynamics of the 
structural transitions for both the first order and continuous cases. Quenches from the 
SQR to the TRN lattice in this system may be accomplished simply by changing d. 
How does the new phase form inside the parent? Is there a possibility of a martensitic 
transition [22l [23]? If so, then of what type? We hope our work stimulates experiments 
designed to answer these questions in the near future. 
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